Mode excitation Monte Carlo simulations of mesoscopically large membranes 
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Solvent-free coarse grained models represent one of the most promising approaches for molecular 
simulations of mesoscopically large membranes. In these models, the size of the simulated membrane 
is limited by the slow relaxation time of longest bending mode. Here, we present a Monte Carlo 
algorithm with update moves in which all the lipids are displaced simultaneously. These collective 
moves result in fast excitation and relaxation of the long wavelength thermal fluctuations. We apply 
the method to simulations of a bilayer membrane of linear size ~ 50 nm and show reduction of the 
relaxation time by two orders of magnitudes when compared to conventional Monte Carlo. 

PACS numbers: 



Biological membranes play a vital role in almost all cel- 
lular phenomena and are fundamental to the organization 
of the cell. Because of their remarkable complexity, com- 
puter models have become essential to the understanding 
of their structure and dynamics. Computer simulations 
of lipid and biological membranes can be broadly clas- 
sified into (i) atomistic models which are limited in the 
size and time of problems they can address by their huge 
computational workload and (ii) coarse grained (CG) 
models that sacrifice most of the atomistic details in or- 
der to explore larger length- and time-scales. The field of 
simplified membrane simulations is more than 20 years 
old and goes back to the work of Kantor et al. on solid 
tethered membranes ^] , which was later extended to sim- 
ulations of fluid membrane by considering dynamically 
triangulated networks 0]. A few years later, molecular 
bead-spring lipid models were developed to elucidate mi- 
celle self-assembly in aqueous environment Q • Recently, 
a new class of CG molecular models have been introduced 
in which bilayer membranes and vesicles are simulated 
without direct representation of an embedding solvent 
0. This is accomplished by constructing intermolecular 
force fields that mimic effects of hydration. The develop- 
ment of implicit-solvent models constitutes an important 
advance in large-scale membrane simulations, considering 
the fact that the number of solvent particles in explicit- 
solvent models is significantly larger than the number of 
lipids. These models now serve as platforms for simula- 
tions of complexes of membranes with proteins [6|, 0| and 
DNA molecules 3- 

Existing implicit solvent CG bilayer models employ an 
extremely simple representation of the lipids as short 
chains consisting of one hydrophilic bead (representing 
the head group) and two hydrophobic beads (represent- 
ing the hydrocarbon tail), connected to each other by 
stiff springs. In earlier works, simulations of membranes 
consisting of ~ 1000 lipids have been presented @. 
Such simulations can be easily performed on a commod- 
ity PC/workstation. A membrane patch of 1000 lipids 
has the linear size of about L ^ 2Q nm (taking the area 
per lipid to be ^ 0.7 nm^), which is at the small-size end 
of the mesoscopic regime. Simulations of larger mem- 
branes would require more memory storage and CPU 
time. The memory needed for simulations of membranes 



containing A^ lO'* lipids is still significantly smaller 
than the memory available on a normal PC. The CPU 
time problem, however, is formidable. For tensionless 
membranes, the relaxation time of the longest bending 
mode scales as r ~ ~ A^^ (see Eq.(|31) below). More- 
over, for particles interacting via short-range forces only, 
the CPU time per Monte Carlo (MC) or Brownian Dy- 
namics time step scales as A^. Therefore, the total CPU 
time of the simulations grows as A''^, from roughly 10 
hours for A^ = 10'^ (on an AMD Opteron 275 processor 
running at 2.2 Ghz) to more than a year (!) for A^ = 10^. 
In this work we propose an improved MC scheme that 
considerably reduces these enormous computing times 
and, thus, permit simulations of membrane-based sys- 
tems on larger length and time scales. 

Recent computer simulations by Reynwar et al. Q 
demonstrate the slow relaxation problem in membrane 
simulations: In this work, the assembly of membrane in- 
clusions by curvature-mediated interactions was studied 
[lo| . Calculating the interaction between a pair of in- 
clusion requires that the equilibrium statistics of ther- 
mal fluctuations on the scale of the object pair sepa- 
ration distance are accurately measured. To access the 
regime of large separations, Reynwar et al. employ a CG 
implicit-solvent model, which permitted them to sim- 
ulate a square membrane of 46080 lipids (the largest 
membrane patch simulated to date) with a linear size 
of about L 130 nm. Unfortunately, the membrane- 
mediated interactions cannot be computed at these large 
spatial scales because the temporal evolution of the cor- 
responding bending modes is extremely slow. There- 
fore, the membrane in this study is decorated with 36 
inclusions initially forming a square lattice with spacing 
(i = L/6 ~ 20 nm, and the calculation of the forces is 
limited to this range. 

The slowing down of single particle update schemes 
arises because the relaxation of large scale fluctuations 
requires a coordinated movement of all the lipids over in- 
creasingly larger distances. In lattice membrane models 
this problem can be solved by using the Fourier represen- 
tation of the membrane height field and updating one, 
randomly chosen, Fourier amplitude at a time [11] . (The 
method was originally proposed for lattice gauge models 
p^.) The efficiency of this method relies the fact that in 
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lattice simulations each Fourier mode of the membrane 
represents a single degree of freedom of the field and, 
therefore, large variations in their amplitudes will not be 
energetically costly and will have reasonable acceptance 
probabilities. Such large scale variations are prohibited 
in off-lattice molecular models by excluded volume inter- 
actions. Nevertheless, MC algorithms with collective up- 
date moves have been recently proposed for simulations 
of simple fluids exhibiting superior performance 

over conventional MC schemes. For molecular simula- 
tions of membranes and interfaces, we consider collective 
moves in which the coordinates of all the lipids are si- 
multaneously updated according to the rule: 



a^old, J/old, Zold 



£i cos{qi 



i=l 



(1) 

where the sum runs over a set of m modes with wavevec- 
tors q = (27r/L) • (711,712); ni,n2 = 0,1,2,.... This 
set includes the modes with the smallest wavenumbers 
n'^ = n1 + ~ Ij 2, 4, 5,8,.. .. The random amplitude 
of the i-th mode in Eq.([T]) is chosen from the interval 
[— £*, -l-e*] (the magnitude of e* is discussed below), while 
ai is a random phase chosen from a uniform distribution 
on [0, 2tt). Because this move is reversed by choosing the 
set of amplitudes {— 6^}, and since the Jacobian of trans- 
formation described by Eq.([T]) is unity, detailed balance 
is satisfied when the proposed mapping Eq.([T]) is com- 
bined with Metropolis acceptance rule: p{o\d new) — 
min(l, exp(— /3Ai<^)), where /3 = {kBT)~^ is the inverse 
temperature and Ai? is the energy difference between the 
"new" and "old" states. 

The height function of the simulated membrane is cal- 
culated by dividing the area into JVP — {L/lY grid cell 
of size / (comparable to the width of the membrane), 
and averaging the height of all the lipids instantaneously 
located within each cell The Fourier transform of 

the discrete height function (defined on the set of points 
{r^}, each of which is located at the center of a grid cell) 
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includes modes, corresponding to ft = 

(ni,n2) ; ni,n2 = (-M/2) -h 1, . . . , M/2 . In con- 
ventional MC simulations, all the modes are equally 
affected by the uncorrelated move attempts. Randomly 
displacing the lipids a vertical distance e within a MC 
time unit, would cause the amplitudes of all the Fourier 
modes @ to change by {SKf ~ M^e^ = {L/lYt^, 
independently of n. At large scales (small n), the 
behaviour of an undulating membrane can be described 
by Helfrich effective surface Hamiltonian which relates 
the elastic energy to the local curvature and the bending 
modulus, K. The power spectrum of the bilayer thermal 
fluctuations [16| 
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strongly depends on n. The conventional MC scheme 
generates diffusive dynamics in Fourier space, where for 
each mode the relaxation time (in MC time units) 



{\hn?) _ kT 



(4) 



The relaxation time of the mode with the smallest 
wavenumber, n — 1, grows as a fourth power of the linear 
size of the system. 

The relaxation of long wavelength modes can be accel- 
erated by introducing collective MC moves which cause 
larger variations in their amplitudes. To eliminate the 
n~'^ dependence of r„ and ensure that all the modes relax 
equally fast, we set the interval from which the random 
amplitudes in Eq.([T]) are chosen to satisfy: e* = A/n^ 
(see Eq.Q). The value of A can be determined empir- 
ically, by employing the usual criterion that the accep- 
tance rate of the moves represented by Eq.((T]) is approxi- 
mately half. Notice, however, that because of the strong 
decrease of e* with n, significant improvement in the re- 
laxation times should be expected only for the longer 
(also the slower) wavelength modes. Therefore, the sum 
in Eq. ^ can be limited to small wavenumbers while the 
relaxation of modes with larger values of n will contin- 
ued to rely on single particle moves. The long wave- 
length modes are efficiently sampled by the new scheme 
because the magnitude of A is independent of the sys- 
tem size. This can be understood by noting that the 
energy cost per unit area of a collective trial move is 
E/L'^ CS'^, where S ~ (A/L) is the induced strain 
and C is the relevant elastic modulus. The total defor- 
mation energy E ^ CA? should be of the order of the 
thermal energy /cbT, yielding A^ ~ {kBT)/C which is 
indeed size- independent. The collective MC moves cause 
the amplitudes of the slow modes included in the sum in 
Eq.© to change by \Shn\ ~ M^A = {L/lf{A/n^) (see 
Eq.©) and, therefore, their relaxation times scale as 



K A2(27r)' 



(5) 



This time does not increase with decreasing the 
wavenumber n and, moreover, grows only ds LP' ^ N 
rather than L"* ~ A'^^. Furthermore, a single collec- 
tive trial move requires the evaluation of 0{N) (short 
range) pair-interactions, which makes them equally CPU 
time as 0{N) single particle trial moves. Therefore, the 
CPU time per MC time unit required in schemes utilizing 
0{N) single-particle and 0(1) collective mode excitation 
trial moves would scale as N'^, which is superior to con- 
ventional MC algorithms whose CPU time grows as N^. 

To demonstrate the validity and efhciency of the new 
algorithm, we carried out simulations using Reynwar et 
al. three-bead lipid model. The details of the intra- and 
intermolecular potentials are given in ref (l7j . In our 
study we set the energy parameter of the Lennard- Jones 
(LJ) potential e = l.OS/csT and the range of the attrac- 
tive tail-tail potential Wc = 1.35cr, where a is the length 




FIG. 1: The normalized distribution functions p of the pro- 
jected area per Upid a/a^, where a is the length parameter of 
the bead-bead LJ potential. 
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FIG. 2; Fluctuation spectrum of a membrane of TV = 1000 
lipids. Results of the conventional and improved algorithms 
are shown by open circles and dashed line, respectively. The 



solid line indicates the asymptotic {\hn\ ) 



power law. 



parameter of the LJ potential. For this choice of the 
parameters, the membrane is in the fluid phase. The 
intermolecular interactions were slightly modified from 
the original model to eliminate the occasional escape of 
lipids from the membrane plane, without affecting the 
rigidity and fluidity of the membrane. To verify that the 
newly proposed mode excitation Monte Carlo (MEMC) 
algorithm works correctly, we used it for MC simulations 
of square membranes with TV = 1000 lipids and com- 
pared the results to those obtained by a conventional 
MC algorithm. The simulations were conducted in the 
constant surface tension ensemble [isj . at vanishing sur- 
face tension. In the conventional algorithm, each MC 
time unit consisted of N displacement move attempts of 
lipids (including changes in the relative coordinates of 
the beads), N rotation move attempts, and two area- 
changing trial moves. The improved MEMC algorithm 
included two additional trial moves per MC time unit in 
which all the modes with wavenumbers rt^ < 8 in Eq.([T]) 
are excited. The (normalized) distribution functions of 
the projected area per lipid, a — 2L^/TV, obtained from 
the conventional and improved simulations are plotted in 
Fig[T] Within negligible computational uncertainties the 
two distribution functions are indistinguishable, which 
confirms that both algorithms generate the same statis- 
tical ensembles. The power spectrum P) of the height 
fluctuations is plotted in Fig. O The conventional and 
improved algorithms give identical results, including the 
asymptotic ^ power law. From Eq.Q (set- 

ting the mesh size to Z = L/8), we calculate the bending 
modulus of the bilayer n ~ SksT, in consistency with 
the values measured in ref [l7j . 

Next, we tested the improvement in computational 
efficiency by simulating larger membranes consisting of 
N — 9000 lipids. The cross sectional area of the simu- 
lation cell was divided into a 16 x 16 grid and the (dis- 
crete) height function was evaluated every 50 MC time 
units. The Fourier transform of height function Q was 
then computed and the amplitudes off all the modes with 



wavenumbers rv^ < 29 were recorded. The relaxation 
times were calculated by fitting the time autocorrelation 
function: C„(At) = (|/i„(i)/i„(t At)|)/(|/i„(t)|2) to a 
double exponential function: C„(At) = a exp{— At / t^) + 
(1 — a) exp(— Af/r^). The double exponential decay has 
been originally conjectured by Seifert and Langer p^ . 
and was recently observed in simulations of Shkulipa et 
al. po| . In our simulations, the dissipation of the bend- 
ing energy accounts for the slow relaxation mechanism 
characterized by t„^i = r„, while the smaller relaxation 
time T„^2 may be associated with intermonolayer friction. 
The latter mechanism was found to play only a relatively 
minor role in the decay of all the investigated modes. 
Due to the large statistical noise and in order to reduce 
the cross correlation between the two relaxation times, 
the fit to a double exponential form was limited to time 
intervals t„_2 <C Ai < t„. The uncertainties in t„ (typ- 
ically ±20 — 25%) were determined by comparing the 
fit results obtained for different fitting intervals. In the 
MEMC algorithm, each MC time unit consisted of (on 
average): TV translations, TV rotations, 2 area-changing, 
and 18 mode (with wavenumbers n? < 13) excitation 
trial moves. The conventional MC algorithm included 
only the first three move types applied with TV : TV : 2 
proportions; however, each MC time unit of the conven- 
tional algorithm consisted of almost 8TV trails in order 
to make the CPU time per MC time unit of both al- 
gorithms the same. The results of our analysis of the 
relaxation times are summarized in Fig. [31 The MEMC 
algorithm eliminates the slowing down of the long wave- 
length modes (solid squares), causing them to relax at 
very similar rates. The relaxation times of the short 
modes which are not excited (open squares) follow the 
T,i ^ n^^ power law (dashed lines), which is also obeyed 
by the modes when the conventional MC algorithm is ap- 
plied (open circles). At small length scales the MEMC 
algorithm is almost 4 times slower than the conventional 
scheme because each MC time unit of the latter includes 



4 











1 1 111. 


le+06 








- 


le+05 










e 

10000 




i 






1000 











J , , I , , III 

1 10 



2 

n 

FIG. 3: Relaxation times of undulatory bending modes as a 
function of the wavenumber n. Conventional MC results are 
plotted by circles (open: results obtained numerically, solid: 
results evaluated by extrapolation). MEMC results are plot- 
ted by squares (open: unexcited modes, solid: excited modes). 
Dashed lines: fits to the r ~ n^* scaling law. 

almost 4 times more single particle moves. The relax- 
ation of the total bilayer area (not shown) , which is quite 
fast, is also slowed down by a factor of 4. The relaxation 
times of the long (excited) modes are considerably re- 
duced and become comparable to the relaxation times of 
the longest modes among those which were not excited by 
the collective update moves. In comparison to the con- 
ventional scheme, the relaxation of the n = 1 modes is 
improved by a factor of about 50, from an estimated one 
year of CPU time to less than a week. The simulations 
extended over a period of about 10 weeks and, therefore, 
our estimates of the long scales slow relaxation times for 
the conventional MC algorithm (solid circles) is based on 
extrapolation of the r„ power law rather than on 

direct numerical evaluation. 

It is interesting to compare the efficiency of the MEMC 
algorithm with alternative computational algorithms for 
constant temperature simulations. MEMC is clearly 
more efficient than constant temperature molecular dy- 
namics (MD) algorithms which at sufficiently large scales 
become Brownian in nature and effectively behave like 
conventional MC simulations [2l|. Improved relaxation 
behavior is achieved when the MD simulations are run 



with a momentum-conserving thermostat |22| that, on 
long length and time scales, reproduce the correct hydro- 
dynamic behavior t [23|. When the CPU time per 
time step is considered, one finds that the computational 
complexity of such simulations grows as iV^ "'^. This is 
better than conventional MC and MD but still inferior to 
the MEMC algorithm whose complexity grows as N'^. In 
the lattice membrane simulations [11,] , the CPU time per 
MC time step grows as iV^ (since there are 0{N) Fourier 
modes and the variation of each is a collective move 
that requires the calculation of 0{N) interaction terms), 
which makes it comparable to MEMC simulations of ten- 
sionless membranes when each lattice point represents a 
microscopic area element of the membrane. However, 
when the membrane is under tension or in the presence 
of an external harmonic potential, the power spectrum 
for small wave numbers are given by P) ~ n^'^ and 
{\hfi p) ~ n°, respectively di3]. Repeating the argument 
that leads from Eq. ^ to Eq. ([5|) , one finds t„ ^ L'^ rather 
than r„ ^ L? in Eq.© and, therefore, the required CPU 
time for MEMC simulations of such membranes would 
grow only linearly with N . This demonstrates that the 
MEMC algorithm is asymptotically faster not only than 
other algorithms for continuum (molecular) membrane 
simulations, but also than the Fourier MC algorithm for 
lattice simulations. 

To summarize, we introduce an improved MC algo- 
rithm for simulations of mesoscopically large membranes. 
The new algorithm utilizes collective update moves that 
lead to fast excitation and relaxation of the long wave- 
length bending modes. The slow relaxation of these 
modes in conventional MC and MD schemes is the most 
severe constraint that limits the size of the simulated 
membranes in solvent-free coarse grained models. The 
efficiency of the new algorithm is demonstrated by sim- 
ulations of a membrane patch of 9000 lipids, where a 50- 
fold decrease in the relaxation time was measured as com- 
pared to a conventional MC algorithm with only single 
particle moves. Implicit solvent bilayer models combined 
with improved sampling techniques, such as the mode 
excitation algorithm presented here, can serve as the ba- 
sis for large scale CG simulations of complexes of bilayer 
membranes with additional biological components. 
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